using System;
using System.Collections.Generic;
using System.Text;
using NUnit.Framework;
using dnAnalytics.Math;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Decomposition;
using dnAnalytics.Exceptions;

namespace dnAnalytics.UnitTests.LinearAlgebra.Decomposition
{
    [TestFixture]
    [Category("Managed")]
    public class SymmetricLevinsonTest
    {
        #region Fields

        // unit testing - order 1

        Vector T1;			// Toeplitz matrix
        Matrix A1;			// Lower triangle matrix
        Vector D1;			// diagonal vector
        Matrix I1;			// inverse matrix
        double Det1;				// exact determinant
        Vector X1;			// RHS vector
        Vector Y1;			// LHS vector
        double Tolerance1;			// allowable tolerance

        // unit testing - order 2

        Vector T2;			// Toeplitz matrix
        Matrix A2;			// Lower triangle matrix
        Vector D2;			// diagonal vector
        Matrix I2;			// inverse matrix
        double Det2;				// exact determinant
        Vector X2;			// RHS vector
        Vector Y2;			// LHS vector
        double Tolerance2;			// allowable tolerance

        // unit testing - order 3

        Vector T3;			// Toeplitz matrix
        Matrix A3;			// Lower triangle matrix
        Vector D3;			// diagonal vector
        Matrix I3;			// inverse matrix
        double Det3;				// exact determinant
        Vector X3;			// RHS vector
        Vector Y3;			// LHS vector
        double Tolerance3;			// allowable tolerance

        // unit testing - order 4

        Vector T4;			// Toeplitz matrix
        Matrix A4;			// Lower triangle matrix
        Vector D4;			// diagonal vector
        Matrix I4;			// inverse matrix
        double Det4;				// exact determinant
        Vector X4;			// RHS vector
        Vector Y4;			// LHS vector
        double Tolerance4;			// allowable tolerance

        // unit testing - order 5

        Vector T5;			// Toeplitz matrix
        Matrix A5;			// Lower triangle matrix
        Vector D5;			// diagonal vector
        Matrix I5;			// inverse matrix
        double Det5;				// exact determinant
        Vector X5;			// RHS vector
        Vector Y5;			// LHS vector
        double Tolerance5;			// allowable tolerance

        // unit testing - order 10

        Vector T10;			// Toeplitz matrix
        Matrix A10;			// Lower triangle matrix
        Vector D10;			// diagonal vector
        Matrix I10;			// inverse matrix
        double Det10;				// exact determinant
        Vector X10;			// RHS vector
        Vector Y10;			// LHS vector
        double Tolerance10;			// allowable tolerance

        #endregion Fields

        #region Test Fixture Setup

        [TestFixtureSetUp]
        public void SetupTestCases()
        {
            // unit testing values - order 1

            T1 = VectorBuilder.CreateVector(1, VectorType.Dense);
            T1[0] = +3.0000000000000000E+000;

            A1 = MatrixBuilder.CreateMatrix(1);
            A1[0, 0] = +1.0000000000000000E+000;

            D1 = VectorBuilder.CreateVector(1, VectorType.Dense);
            D1[0] = +3.3333333333333331E-001;

            Det1 = +3.0000000000000000E+000;

            I1 = MatrixBuilder.CreateMatrix(1);
            I1[0, 0] = +3.3333333333333331E-001;

            X1 = VectorBuilder.CreateVector(1, VectorType.Dense);
            X1[0] = +1.0000000000000000E+000;

            Y1 = VectorBuilder.CreateVector(1, VectorType.Dense);
            Y1[0] = +3.0000000000000000E+000;

            Tolerance1 = 1.0E-015;

            // unit testing values - order 2

            T2 = VectorBuilder.CreateVector(2, VectorType.Dense);
            T2[0] = +1.0000000000000000E+001;
            T2[1] = +4.5000000000000000E+000;

            A2 = MatrixBuilder.CreateMatrix(2, MatrixType.Dense);
            A2[0, 0] = +1.0000000000000000E+000;
            A2[1, 0] = -4.5000000000000001E-001;
            A2[1, 1] = +1.0000000000000000E+000;

            D2 = VectorBuilder.CreateVector(2, VectorType.Dense);
            D2[0] = +1.0000000000000001E-001;
            D2[1] = +1.2539184952978055E-001;

            Det2 = +7.9750000000000000E+001;

            I2 = MatrixBuilder.CreateMatrix(2, MatrixType.Dense);
            I2[0, 0] = +1.2539184952978055E-001;
            I2[0, 1] = -5.6426332288401257E-002;
            I2[1, 0] = -5.6426332288401257E-002;
            I2[1, 1] = +1.2539184952978055E-001;

            X2 = VectorBuilder.CreateVector(2, VectorType.Dense);
            X2[0] = +1.0000000000000000E+000;
            X2[1] = +2.0000000000000000E+000;

            Y2 = VectorBuilder.CreateVector(2, VectorType.Dense);
            Y2[0] = +1.9000000000000000E+001;
            Y2[1] = +2.4500000000000000E+001;

            Tolerance2 = 1.0E-015;

            // unit testing values - order 3

            T3 = VectorBuilder.CreateVector(3, VectorType.Dense);
            T3[0] = +4.0000000000000000E+000;
            T3[1] = +1.5000000000000000E+000;
            T3[2] = +6.6666666666666663E-001;

            A3 = MatrixBuilder.CreateMatrix(3, MatrixType.Dense);
            A3[0, 0] = +1.0000000000000000E+000;
            A3[1, 0] = -3.7500000000000000E-001;
            A3[1, 1] = +1.0000000000000000E+000;
            A3[2, 0] = -3.0303030303030304E-002;
            A3[2, 1] = -3.6363636363636365E-001;
            A3[2, 2] = +1.0000000000000000E+000;

            D3 = VectorBuilder.CreateVector(3, VectorType.Dense);
            D3[0] = +2.5000000000000000E-001;
            D3[1] = +2.9090909090909089E-001;
            D3[2] = +2.9117647058823531E-001;

            Det3 = +4.7222222222222221E+001;

            I3 = MatrixBuilder.CreateMatrix(3, MatrixType.Dense);
            I3[0, 0] = +2.9117647058823531E-001;
            I3[0, 1] = -1.0588235294117647E-001;
            I3[0, 2] = -8.8235294117647058E-003;
            I3[1, 0] = -1.0588235294117647E-001;
            I3[1, 1] = +3.2941176470588235E-001;
            I3[1, 2] = -1.0588235294117647E-001;
            I3[2, 0] = -8.8235294117647058E-003;
            I3[2, 1] = -1.0588235294117647E-001;
            I3[2, 2] = +2.9117647058823531E-001;

            X3 = VectorBuilder.CreateVector(3, VectorType.Dense);
            X3[0] = +1.0000000000000000E+000;
            X3[1] = +2.0000000000000000E+000;
            X3[2] = +3.0000000000000000E+000;

            Y3 = VectorBuilder.CreateVector(3, VectorType.Dense);
            Y3[0] = +9.0000000000000000E+000;
            Y3[1] = +1.4000000000000000E+001;
            Y3[2] = +1.5666666666666666E+001;

            Tolerance3 = 3.0E-015;

            // unit testing values - order 4

            T4 = VectorBuilder.CreateVector(4, VectorType.Dense);
            T4[0] = +4.0000000000000000E+000;
            T4[1] = +1.5000000000000000E+000;
            T4[2] = +6.6666666666666663E-001;
            T4[3] = +2.5000000000000000E-001;

            A4 = MatrixBuilder.CreateMatrix(4, MatrixType.Dense);
            A4[0, 0] = +1.0000000000000000E+000;
            A4[1, 0] = -3.7500000000000000E-001;
            A4[1, 1] = +1.0000000000000000E+000;
            A4[2, 0] = -3.0303030303030304E-002;
            A4[2, 1] = -3.6363636363636365E-001;
            A4[2, 2] = +1.0000000000000000E+000;
            A4[3, 0] = +1.1029411764705883E-002;
            A4[3, 1] = -3.4313725490196081E-002;
            A4[3, 2] = -3.6397058823529410E-001;
            A4[3, 3] = +1.0000000000000000E+000;

            D4 = VectorBuilder.CreateVector(4, VectorType.Dense);
            D4[0] = +2.5000000000000000E-001;
            D4[1] = +2.9090909090909089E-001;
            D4[2] = +2.9117647058823531E-001;
            D4[3] = +2.9121189591078067E-001;

            Det4 = +1.6215760030864197E+002;

            I4 = MatrixBuilder.CreateMatrix(4, MatrixType.Dense);
            I4[0, 0] = +2.9121189591078067E-001;
            I4[0, 1] = -1.0599256505576209E-001;
            I4[0, 2] = -9.9925650557620826E-003;
            I4[0, 3] = +3.2118959107806690E-003;
            I4[1, 0] = -1.0599256505576209E-001;
            I4[1, 1] = +3.2975464684014871E-001;
            I4[1, 2] = -1.0224535315985130E-001;
            I4[1, 3] = -9.9925650557620826E-003;
            I4[2, 0] = -9.9925650557620826E-003;
            I4[2, 1] = -1.0224535315985130E-001;
            I4[2, 2] = +3.2975464684014871E-001;
            I4[2, 3] = -1.0599256505576209E-001;
            I4[3, 0] = +3.2118959107806690E-003;
            I4[3, 1] = -9.9925650557620826E-003;
            I4[3, 2] = -1.0599256505576209E-001;
            I4[3, 3] = +2.9121189591078067E-001;

            X4 = VectorBuilder.CreateVector(4, VectorType.Dense);
            X4[0] = +1.0000000000000000E+000;
            X4[1] = +2.0000000000000000E+000;
            X4[2] = +3.0000000000000000E+000;
            X4[3] = +4.0000000000000000E+000;

            Y4 = VectorBuilder.CreateVector(4, VectorType.Dense);
            Y4[0] = +1.0000000000000000E+001;
            Y4[1] = +1.6666666666666668E+001;
            Y4[2] = +2.1666666666666668E+001;
            Y4[3] = +2.2083333333333332E+001;

            Tolerance4 = 4.0E-015;

            // unit testing values - order 5

            T5 = VectorBuilder.CreateVector(5, VectorType.Dense);
            T5[0] = +1.0000000000000000E+000;
            T5[1] = +5.0000000000000000E-001;
            T5[2] = +3.3333333333333331E-001;
            T5[3] = +2.5000000000000000E-001;
            T5[4] = +2.0000000000000001E-001;

            A5 = MatrixBuilder.CreateMatrix(5, MatrixType.Dense);
            A5[0, 0] = +1.0000000000000000E+000;
            A5[1, 0] = -5.0000000000000000E-001;
            A5[1, 1] = +1.0000000000000000E+000;
            A5[2, 0] = -1.1111111111111110E-001;
            A5[2, 1] = -4.4444444444444442E-001;
            A5[2, 2] = +1.0000000000000000E+000;
            A5[3, 0] = -6.2500000000000000E-002;
            A5[3, 1] = -8.3333333333333329E-002;
            A5[3, 2] = -4.3750000000000000E-001;
            A5[3, 3] = +1.0000000000000000E+000;
            A5[4, 0] = -4.2823529411764705E-002;
            A5[4, 1] = -4.3764705882352942E-002;
            A5[4, 2] = -7.9764705882352946E-002;
            A5[4, 3] = -4.3482352941176472E-001;
            A5[4, 4] = +1.0000000000000000E+000;

            D5 = VectorBuilder.CreateVector(5, VectorType.Dense);
            D5[0] = +1.0000000000000000E+000;
            D5[1] = +1.3333333333333333E+000;
            D5[2] = +1.3500000000000001E+000;
            D5[3] = +1.3552941176470588E+000;
            D5[4] = +1.3577840963547489E+000;

            Det5 = +3.0190007716049383E-001;

            I5 = MatrixBuilder.CreateMatrix(5, MatrixType.Dense);
            I5[0, 0] = +1.3577840963547489E+000;
            I5[0, 1] = -5.9039647295613562E-001;
            I5[0, 2] = -1.0830324909747292E-001;
            I5[0, 3] = -5.9423021628701958E-002;
            I5[0, 4] = -5.8145107185073958E-002;
            I5[1, 0] = -5.9039647295613562E-001;
            I5[1, 1] = +1.6120123957701031E+000;
            I5[1, 2] = -5.4584837545126352E-001;
            I5[1, 3] = -8.7102648477684425E-002;
            I5[1, 4] = -5.9423021628701958E-002;
            I5[2, 0] = -1.0830324909747292E-001;
            I5[2, 1] = -5.4584837545126352E-001;
            I5[2, 2] = +1.6180505415162454E+000;
            I5[2, 3] = -5.4584837545126352E-001;
            I5[2, 4] = -1.0830324909747292E-001;
            I5[3, 0] = -5.9423021628701958E-002;
            I5[3, 1] = -8.7102648477684425E-002;
            I5[3, 2] = -5.4584837545126352E-001;
            I5[3, 3] = +1.6120123957701031E+000;
            I5[3, 4] = -5.9039647295613562E-001;
            I5[4, 0] = -5.8145107185073958E-002;
            I5[4, 1] = -5.9423021628701958E-002;
            I5[4, 2] = -1.0830324909747292E-001;
            I5[4, 3] = -5.9039647295613562E-001;
            I5[4, 4] = +1.3577840963547489E+000;

            X5 = VectorBuilder.CreateVector(5, VectorType.Dense);
            X5[0] = +1.0000000000000000E+000;
            X5[1] = +2.0000000000000000E+000;
            X5[2] = +3.0000000000000000E+000;
            X5[3] = +4.0000000000000000E+000;
            X5[4] = +5.0000000000000000E+000;

            Y5 = VectorBuilder.CreateVector(5, VectorType.Dense);
            Y5[0] = +5.0000000000000000E+000;
            Y5[1] = +6.5833333333333330E+000;
            Y5[2] = +8.0000000000000000E+000;
            Y5[3] = +8.9166666666666661E+000;
            Y5[4] = +8.6999999999999993E+000;

            Tolerance5 = 5.0E-015;

            // unit testing values - order 10

            T10 = VectorBuilder.CreateVector(10, VectorType.Dense);
            T10[0] = +1.0000000000000000E+001;
            T10[1] = +5.0000000000000000E+000;
            T10[2] = +3.3333333333333335E+000;
            T10[3] = +2.5000000000000000E+000;
            T10[4] = +2.0000000000000000E+000;
            T10[5] = +1.6666666666666667E+000;
            T10[6] = +1.4285714285714286E+000;
            T10[7] = +1.2500000000000000E+000;
            T10[8] = +1.1111111111111112E+000;
            T10[9] = +1.0000000000000000E+000;

            A10 = MatrixBuilder.CreateMatrix(10, MatrixType.Dense);
            A10[0, 0] = +1.0000000000000000E+000;
            A10[1, 0] = -5.0000000000000000E-001;
            A10[1, 1] = +1.0000000000000000E+000;
            A10[2, 0] = -1.1111111111111110E-001;
            A10[2, 1] = -4.4444444444444442E-001;
            A10[2, 2] = +1.0000000000000000E+000;
            A10[3, 0] = -6.2500000000000000E-002;
            A10[3, 1] = -8.3333333333333329E-002;
            A10[3, 2] = -4.3750000000000000E-001;
            A10[3, 3] = +1.0000000000000000E+000;
            A10[4, 0] = -4.2823529411764705E-002;
            A10[4, 1] = -4.3764705882352942E-002;
            A10[4, 2] = -7.9764705882352946E-002;
            A10[4, 3] = -4.3482352941176472E-001;
            A10[4, 4] = +1.0000000000000000E+000;
            A10[5, 0] = -3.2262015058091863E-002;
            A10[5, 1] = -2.8795246158269703E-002;
            A10[5, 2] = -4.1191335740072200E-002;
            A10[5, 3] = -7.8352768282163504E-002;
            A10[5, 4] = -4.3344195606104169E-001;
            A10[5, 5] = +1.0000000000000000E+000;
            A10[6, 0] = -2.5714690745207582E-002;
            A10[6, 1] = -2.1116189201984319E-002;
            A10[6, 2] = -2.6780428952862961E-002;
            A10[6, 3] = -4.0132113280134232E-002;
            A10[6, 4] = -7.7612307432271474E-002;
            A10[6, 5] = -4.3261234832100565E-001;
            A10[6, 6] = +1.0000000000000000E+000;
            A10[7, 0] = -2.1279070818048949E-002;
            A10[7, 1] = -1.6509101948522443E-002;
            A10[7, 2] = -1.9464671415780829E-002;
            A10[7, 3] = -2.5926454872297021E-002;
            A10[7, 4] = -3.9562250635908527E-002;
            A10[7, 5] = -7.7162974546835128E-002;
            A10[7, 6] = -4.3206516359557418E-001;
            A10[7, 7] = +1.0000000000000000E+000;
            A10[8, 0] = -1.8086542121485090E-002;
            A10[8, 1] = -1.3464506037451252E-002;
            A10[8, 2] = -1.5113490559162029E-002;
            A10[8, 3] = -1.8749127103233718E-002;
            A10[8, 4] = -2.5457534954188437E-002;
            A10[8, 5] = -3.9210202036466142E-002;
            A10[8, 6] = -7.6864381979055291E-002;
            A10[8, 7] = -4.3168029878491743E-001;
            A10[8, 8] = +1.0000000000000000E+000;
            A10[9, 0] = -1.5685263462054830E-002;
            A10[9, 1] = -1.1315522903665111E-002;
            A10[9, 2] = -1.2258867955261750E-002;
            A10[9, 3] = -1.4498468209819659E-002;
            A10[9, 4] = -1.8349818960382801E-002;
            A10[9, 5] = -2.5163449955890666E-002;
            A10[9, 6] = -3.8973142955214410E-002;
            A10[9, 7] = -7.6653187654471433E-002;
            A10[9, 8] = -4.3139660660662438E-001;
            A10[9, 9] = +1.0000000000000000E+000;

            D10 = VectorBuilder.CreateVector(10, VectorType.Dense);
            D10[0] = +1.0000000000000001E-001;
            D10[1] = +1.3333333333333333E-001;
            D10[2] = +1.3500000000000001E-001;
            D10[3] = +1.3552941176470587E-001;
            D10[4] = +1.3577840963547491E-001;
            D10[5] = +1.3591988015945386E-001;
            D10[6] = +1.3600981601378359E-001;
            D10[7] = +1.3607142900101707E-001;
            D10[8] = +1.3611595566158571E-001;
            D10[9] = +1.3614945216955071E-001;

            Det10 = +6.4761680151808619E+008;

            I10 = MatrixBuilder.CreateMatrix(10, MatrixType.Dense);
            I10[0, 0] = +1.3614945216955071E-001;
            I10[0, 1] = -5.8734411657295094E-002;
            I10[0, 2] = -1.0436289506206054E-002;
            I10[0, 3] = -5.3061720626780265E-003;
            I10[0, 4] = -3.4259899261904190E-003;
            I10[0, 5] = -2.4983177988665531E-003;
            I10[0, 6] = -1.9739585040645931E-003;
            I10[0, 7] = -1.6690381563277477E-003;
            I10[0, 8] = -1.5406022443460088E-003;
            I10[0, 9] = -2.1355400274938358E-003;
            I10[1, 0] = -5.8734411657295094E-002;
            I10[1, 1] = +1.6145378154157938E-001;
            I10[1, 2] = -5.4256396530846282E-002;
            I10[1, 3] = -8.1734041875061048E-003;
            I10[1, 4] = -3.8591736934504169E-003;
            I10[1, 5] = -2.3874108784216197E-003;
            I10[1, 6] = -1.7004963532414416E-003;
            I10[1, 7] = -1.3371698139059260E-003;
            I10[1, 8] = -1.1681235264574500E-003;
            I10[1, 9] = -1.5406022443460088E-003;
            I10[2, 0] = -1.0436289506206054E-002;
            I10[2, 1] = -5.4256396530846282E-002;
            I10[2, 2] = +1.6223632367953364E-001;
            I10[2, 3] = -5.3868547567483924E-002;
            I10[2, 4] = -7.9331275114551281E-003;
            I10[2, 5] = -3.6959394426671018E-003;
            I10[2, 6] = -2.2748675342649503E-003;
            I10[2, 7] = -1.6326013697479989E-003;
            I10[2, 8] = -1.3371698139059260E-003;
            I10[2, 9] = -1.6690381563277477E-003;
            I10[3, 0] = -5.3061720626780265E-003;
            I10[3, 1] = -8.1734041875061048E-003;
            I10[3, 2] = -5.3868547567483924E-002;
            I10[3, 3] = +1.6242266136350664E-001;
            I10[3, 4] = -5.3759224468977879E-002;
            I10[3, 5] = -7.8663867627389306E-003;
            I10[3, 6] = -3.6610068338217565E-003;
            I10[3, 7] = -2.2748675342649503E-003;
            I10[3, 8] = -1.7004963532414416E-003;
            I10[3, 9] = -1.9739585040645931E-003;
            I10[4, 0] = -3.4259899261904190E-003;
            I10[4, 1] = -3.8591736934504169E-003;
            I10[4, 2] = -7.9331275114551281E-003;
            I10[4, 3] = -5.3759224468977879E-002;
            I10[4, 4] = +1.6248025171494504E-001;
            I10[4, 5] = -5.3732579955257077E-002;
            I10[4, 6] = -7.8663867627389306E-003;
            I10[4, 7] = -3.6959394426671018E-003;
            I10[4, 8] = -2.3874108784216197E-003;
            I10[4, 9] = -2.4983177988665531E-003;
            I10[5, 0] = -2.4983177988665531E-003;
            I10[5, 1] = -2.3874108784216197E-003;
            I10[5, 2] = -3.6959394426671018E-003;
            I10[5, 3] = -7.8663867627389306E-003;
            I10[5, 4] = -5.3732579955257077E-002;
            I10[5, 5] = +1.6248025171494504E-001;
            I10[5, 6] = -5.3759224468977879E-002;
            I10[5, 7] = -7.9331275114551281E-003;
            I10[5, 8] = -3.8591736934504169E-003;
            I10[5, 9] = -3.4259899261904190E-003;
            I10[6, 0] = -1.9739585040645931E-003;
            I10[6, 1] = -1.7004963532414416E-003;
            I10[6, 2] = -2.2748675342649503E-003;
            I10[6, 3] = -3.6610068338217565E-003;
            I10[6, 4] = -7.8663867627389306E-003;
            I10[6, 5] = -5.3759224468977879E-002;
            I10[6, 6] = +1.6242266136350664E-001;
            I10[6, 7] = -5.3868547567483924E-002;
            I10[6, 8] = -8.1734041875061048E-003;
            I10[6, 9] = -5.3061720626780265E-003;
            I10[7, 0] = -1.6690381563277477E-003;
            I10[7, 1] = -1.3371698139059260E-003;
            I10[7, 2] = -1.6326013697479989E-003;
            I10[7, 3] = -2.2748675342649503E-003;
            I10[7, 4] = -3.6959394426671018E-003;
            I10[7, 5] = -7.9331275114551281E-003;
            I10[7, 6] = -5.3868547567483924E-002;
            I10[7, 7] = +1.6223632367953364E-001;
            I10[7, 8] = -5.4256396530846282E-002;
            I10[7, 9] = -1.0436289506206054E-002;
            I10[8, 0] = -1.5406022443460088E-003;
            I10[8, 1] = -1.1681235264574500E-003;
            I10[8, 2] = -1.3371698139059260E-003;
            I10[8, 3] = -1.7004963532414416E-003;
            I10[8, 4] = -2.3874108784216197E-003;
            I10[8, 5] = -3.8591736934504169E-003;
            I10[8, 6] = -8.1734041875061048E-003;
            I10[8, 7] = -5.4256396530846282E-002;
            I10[8, 8] = +1.6145378154157938E-001;
            I10[8, 9] = -5.8734411657295094E-002;
            I10[9, 0] = -2.1355400274938358E-003;
            I10[9, 1] = -1.5406022443460088E-003;
            I10[9, 2] = -1.6690381563277477E-003;
            I10[9, 3] = -1.9739585040645931E-003;
            I10[9, 4] = -2.4983177988665531E-003;
            I10[9, 5] = -3.4259899261904190E-003;
            I10[9, 6] = -5.3061720626780265E-003;
            I10[9, 7] = -1.0436289506206054E-002;
            I10[9, 8] = -5.8734411657295094E-002;
            I10[9, 9] = +1.3614945216955071E-001;

            X10 = VectorBuilder.CreateVector(10, VectorType.Dense);
            X10[0] = +1.0000000000000000E+000;
            X10[1] = +2.0000000000000000E+000;
            X10[2] = +3.0000000000000000E+000;
            X10[3] = +4.0000000000000000E+000;
            X10[4] = +5.0000000000000000E+000;
            X10[5] = +6.0000000000000000E+000;
            X10[6] = +7.0000000000000000E+000;
            X10[7] = +8.0000000000000000E+000;
            X10[8] = +9.0000000000000000E+000;
            X10[9] = +1.0000000000000000E+001;

            Y10 = VectorBuilder.CreateVector(10, VectorType.Dense);
            Y10[0] = +1.0000000000000000E+002;
            Y10[1] = +1.2328968253968254E+002;
            Y10[2] = +1.4769047619047620E+002;
            Y10[3] = +1.7195238095238096E+002;
            Y10[4] = +1.9500000000000000E+002;
            Y10[5] = +2.1566666666666666E+002;
            Y10[6] = +2.3242857142857142E+002;
            Y10[7] = +2.4294047619047620E+002;
            Y10[8] = +2.4289682539682539E+002;
            Y10[9] = +2.2218650793650792E+002;

            Tolerance10 = 2.0E-014;
        }

        #endregion Test Fixture Setup

        #region Null Parameter Tests for Constructor

        // Test constructor with a null parameter
        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void NullParameterTestforConstructor()
        {
            double[] data = null;
            SymmetricLevinson dsl = new SymmetricLevinson(data);
        }

        // Test constructor with a null parameter
        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void NullParameterTestforConstructor2()
        {
            Vector data = null;
            SymmetricLevinson dsl = new SymmetricLevinson(data);
        }

        #endregion  Null Parameter Tests for Constructor

        #region Zero Length Vector Tests for Constructor

        [Test]
        [ExpectedException(typeof(System.RankException))]
        public void ZeroLengthVectorTestsforConstructor()
        {
            double[] dv = new double[0];
            SymmetricLevinson dsl = new SymmetricLevinson(dv);
        }

        #endregion Zero Length Vector Tests for Constructor

        #region ToeplitzMatrix Member Test

        // check get matrix
        [Test]
        public void ToeplitzMatrixMember()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Matrix dsldm = dsl.ToeplitzMatrix();

            for (int row = 0; row < T5.Count; row++)
            {
                for (int column = 0; column < T5.Count; column++)
                {
                    if (column < row)
                    {
                        Assert.IsTrue(dsldm[row, column] == T5[row - column]);
                    }
                    else
                    {
                        Assert.AreEqual(dsldm[row, column], T5[column - row], Constants.ACCEPTABLE_ERROR);
                    }
                }
            }
        }

        #endregion ToeplitzMatrix Member Test

        #region ToeplitzMatrix Member Result Test

        // check get matrix
        [Test]
        public void ToeplitzMatrixMemberResult()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Matrix dsldm = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.ToeplitzMatrix(dsldm);

            for (int row = 0; row < T5.Count; row++)
            {
                for (int column = 0; column < T5.Count; column++)
                {
                    if (column < row)
                    {
                        Assert.IsTrue(dsldm[row, column] == T5[row - column]);
                    }
                    else
                    {
                        Assert.AreEqual(dsldm[row, column], T5[column - row], Constants.ACCEPTABLE_ERROR);
                    }
                }
            }
        }

        #endregion ToeplitzMatrix Member Test

        #region ToeplitzMatrix Member Result Error Tests

        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void ToeplitzMatrixMemberResultNull()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Matrix dsldm = null;
            dsl.ToeplitzMatrix(dsldm);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void ToeplitzMatrixMemberResultNotConformable()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Matrix dsldm = MatrixBuilder.CreateMatrix(dsl.Order + 1);
            dsl.ToeplitzMatrix(dsldm);
        }

        #endregion ToeplitzMatrix Member Result Error Tests

        #region Upper Result Error Tests

        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void UpperResultNull()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Matrix upper = null;
            dsl.Upper(upper);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void UpperResultNotConformable()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Matrix upper = MatrixBuilder.CreateMatrix(dsl.Order + 1);
            dsl.Upper(upper);
        }

        #endregion Upper Result Error Tests

        #region Lower Result Error Tests

        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void LowerResultNull()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Matrix lower = null;
            dsl.Lower(lower);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void LowerResultNotConformable()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Matrix lower = MatrixBuilder.CreateMatrix(dsl.Order + 1);
            dsl.Lower(lower);
        }

        #endregion Lower Result Error Tests

        #region Diagonal Result Error Tests

        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void DiagonalResultNull()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Matrix diagonal = null;
            dsl.Diagonal(diagonal);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void DiagonalNotConformable()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Matrix diagonal = MatrixBuilder.CreateMatrix(dsl.Order + 1);
            dsl.Diagonal(diagonal);
        }

        #endregion Diagonal Result Error Tests

        #region Order Property Test

        // test order property
        [Test]
        public void OrderProperty()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Assert.IsTrue(dsl.Order == 5);
        }

        #endregion Order Property Test

        #region Decomposition Test 1

        // test the UDL factorisation for case 1
        [Test]
        public void Decomposition1()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T1);
            Matrix U = dsl.Upper();
            Matrix D = dsl.Diagonal();
            Matrix L = dsl.Lower();

            // check U is the transpose of L
            Matrix LT = L.Transpose();
            for (i = 0; i < LT.Rows; i++)
            {
                for (j = 0; j < LT.Columns; j++)
                {
                    Assert.AreEqual(U[i, j], LT[i, j]);
                }
            }

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    e = System.Math.Abs((A1[i, j] - L[i, j]) / A1[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((D1[i] - D[i, i]) / D1[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 1

        #region Decomposition Test 2

        // test the UDL factorisation for case 2
        [Test]
        public void Decomposition2()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T2);
            Matrix U = dsl.Upper();
            Matrix D = dsl.Diagonal();
            Matrix L = dsl.Lower();

            // check U is the transpose of L
            Matrix LT = L.Transpose();
            for (i = 0; i < LT.Rows; i++)
            {
                for (j = 0; j < LT.Columns; j++)
                {
                    Assert.AreEqual(U[i, j], LT[i, j]);
                }
            }

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    e = System.Math.Abs((A2[i, j] - L[i, j]) / A2[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((D2[i] - D[i, i]) / D2[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 2

        #region Decomposition Test 3

        // test the UDL factorisation for case 3
        [Test]
        public void Decomposition3()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T3);
            Matrix U = dsl.Upper();
            Matrix D = dsl.Diagonal();
            Matrix L = dsl.Lower();

            // check U is the transpose of L
            Matrix LT = L.Transpose();
            for (i = 0; i < LT.Rows; i++)
            {
                for (j = 0; j < LT.Columns; j++)
                {
                    Assert.AreEqual(U[i, j], LT[i, j]);
                }
            }

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    e = System.Math.Abs((A3[i, j] - L[i, j]) / A3[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((D3[i] - D[i, i]) / D3[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 3

        #region Decomposition Test 4

        // test the UDL factorisation for case 4
        [Test]
        public void Decomposition4()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T4);
            Matrix U = dsl.Upper();
            Matrix D = dsl.Diagonal();
            Matrix L = dsl.Lower();

            // check U is the transpose of L
            Matrix LT = L.Transpose();
            for (i = 0; i < LT.Rows; i++)
            {
                for (j = 0; j < LT.Columns; j++)
                {
                    Assert.AreEqual(U[i, j], LT[i, j]);
                }
            }

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    e = System.Math.Abs((A4[i, j] - L[i, j]) / A4[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((D4[i] - D[i, i]) / D4[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 4

        #region Decomposition Test 5

        // test the UDL factorisation for case 5
        [Test]
        public void Decomposition5()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Matrix U = dsl.Upper();
            Matrix D = dsl.Diagonal();
            Matrix L = dsl.Lower();

            // check U is the transpose of L
            Matrix LT = L.Transpose();
            for (i = 0; i < LT.Rows; i++)
            {
                for (j = 0; j < LT.Columns; j++)
                {
                    Assert.AreEqual(U[i, j], LT[i, j]);
                }
            }

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    e = System.Math.Abs((A5[i, j] - L[i, j]) / A5[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((D5[i] - D[i, i]) / D5[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 5

        #region Decomposition Test 10

        // test the UDL factorisation for case 10
        [Test]
        public void Decomposition10()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            Matrix U = dsl.Upper();
            Matrix D = dsl.Diagonal();
            Matrix L = dsl.Lower();

            // check U is the transpose of L
            Matrix LT = L.Transpose();
            for (i = 0; i < LT.Rows; i++)
            {
                for (j = 0; j < LT.Columns; j++)
                {
                    Assert.AreEqual(U[i, j], LT[i, j]);
                }
            }

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    e = System.Math.Abs((A10[i, j] - L[i, j]) / A10[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((D10[i] - D[i, i]) / D10[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 10

        #region Decomposition Result Test 1

        // test the UDL factorisation for case 1
        [Test]
        public void DecompositionResult1()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T1);
            Matrix U = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Upper(U);
            Matrix D = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Diagonal(D);
            Matrix L = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Lower(L);

            // check U is the transpose of L
            Matrix LT = L.Transpose();
            for (i = 0; i < LT.Rows; i++)
            {
                for (j = 0; j < LT.Columns; j++)
                {
                    Assert.AreEqual(U[i, j], LT[i, j]);
                }
            }

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    e = System.Math.Abs((A1[i, j] - L[i, j]) / A1[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((D1[i] - D[i, i]) / D1[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Result Test 1

        #region Decomposition Result Test 2

        // test the UDL factorisation for case 2
        [Test]
        public void DecompositionResult2()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T2);
            Matrix U = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Upper(U);
            Matrix D = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Diagonal(D);
            Matrix L = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Lower(L);

            // check U is the transpose of L
            Matrix LT = L.Transpose();
            for (i = 0; i < LT.Rows; i++)
            {
                for (j = 0; j < LT.Columns; j++)
                {
                    Assert.AreEqual(U[i, j], LT[i, j]);
                }
            }

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    e = System.Math.Abs((A2[i, j] - L[i, j]) / A2[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((D2[i] - D[i, i]) / D2[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Result Test 2

        #region Decomposition Result Test 3

        // test the UDL factorisation for case 3
        [Test]
        public void DecompositionResult3()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T3);
            Matrix U = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Upper(U);
            Matrix D = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Diagonal(D);
            Matrix L = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Lower(L);

            // check U is the transpose of L
            Matrix LT = L.Transpose();
            for (i = 0; i < LT.Rows; i++)
            {
                for (j = 0; j < LT.Columns; j++)
                {
                    Assert.AreEqual(U[i, j], LT[i, j]);
                }
            }

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    e = System.Math.Abs((A3[i, j] - L[i, j]) / A3[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((D3[i] - D[i, i]) / D3[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 3

        #region Decomposition Result Test 4

        // test the UDL factorisation for case 4
        [Test]
        public void DecompositionResult4()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T4);
            Matrix U = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Upper(U);
            Matrix D = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Diagonal(D);
            Matrix L = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Lower(L); ;

            // check U is the transpose of L
            Matrix LT = L.Transpose();
            for (i = 0; i < LT.Rows; i++)
            {
                for (j = 0; j < LT.Columns; j++)
                {
                    Assert.AreEqual(U[i, j], LT[i, j]);
                }
            }

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    e = System.Math.Abs((A4[i, j] - L[i, j]) / A4[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((D4[i] - D[i, i]) / D4[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Result Test 4

        #region Decomposition Result Test 5

        // test the UDL factorisation for case 5
        [Test]
        public void DecompositionResult5()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Matrix U = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Upper(U);
            Matrix D = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Diagonal(D);
            Matrix L = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Lower(L);

            // check U is the transpose of L
            Matrix LT = L.Transpose();
            for (i = 0; i < LT.Rows; i++)
            {
                for (j = 0; j < LT.Columns; j++)
                {
                    Assert.AreEqual(U[i, j], LT[i, j]);
                }
            }

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    e = System.Math.Abs((A5[i, j] - L[i, j]) / A5[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((D5[i] - D[i, i]) / D5[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Result Test 5

        #region Decomposition Result Test 10

        // test the UDL factorisation for case 10
        [Test]
        public void DecompositionResult10()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            Matrix U = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Upper(U);
            Matrix D = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Diagonal(D);
            Matrix L = MatrixBuilder.CreateMatrix(dsl.Order, MatrixType.Dense);
            dsl.Lower(L);

            // check U is the transpose of L
            Matrix LT = L.Transpose();
            for (i = 0; i < LT.Rows; i++)
            {
                for (j = 0; j < LT.Columns; j++)
                {
                    Assert.AreEqual(U[i, j], LT[i, j]);
                }
            }

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    e = System.Math.Abs((A10[i, j] - L[i, j]) / A10[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((D10[i] - D[i, i]) / D10[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Result Test 10

        #region Positve Definite Test 1

        // check the matrix is not positve definite
        [Test]
        public void PositveDefiniteProperty1()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T1);
            Assert.IsTrue(dsl.IsPositiveDefinite());
        }

        #endregion Positve Definite Test 1

        #region Positve Definite Test 2

        // check the matrix is not positve definite
        [Test]
        public void PositveDefiniteProperty2()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T2);
            Assert.IsTrue(dsl.IsPositiveDefinite());
        }

        #endregion Positve Definite Test 2

        #region Positve Definite Test 3

        // check the matrix is not positve definite
        [Test]
        public void PositveDefiniteProperty3()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T3);
            Assert.IsTrue(dsl.IsPositiveDefinite());
        }

        #endregion Positve Definite Test 3

        #region Positve Definite Test 4

        // check the matrix is not positve definite
        [Test]
        public void PositveDefiniteProperty4()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T4);
            Assert.IsTrue(dsl.IsPositiveDefinite());
        }

        #endregion Positve Definite Test 4

        #region Positve Definite Test 5

        // check the matrix is not positve definite
        [Test]
        public void PositveDefiniteProperty5()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Assert.IsTrue(dsl.IsPositiveDefinite());
        }

        #endregion Positve Definite Test 5

        #region Positve Definite Test 10

        // check the matrix is not positve definite
        [Test]
        public void PositveDefiniteProperty10()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            Assert.IsTrue(dsl.IsPositiveDefinite());
        }

        #endregion Positve Definite Test 10

        #region Positve Definite Test - Negative Case

        // check the matrix is not positve definite
        [Test]
        public void PositveDefiniteProperty()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(new double[] { 1.0, 2.0, 3.0, 4.0, 5.0 });
            Assert.IsFalse(dsl.IsPositiveDefinite());
        }

        #endregion Positve Definite Test - Negative Case

        #region Singularity Test 1

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty1()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T1);
            Assert.IsFalse(dsl.IsSingular());
        }

        #endregion Singularity Test 1

        #region Singularity Test 2

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty2()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T2);
            Assert.IsFalse(dsl.IsSingular());
        }

        #endregion Singularity Test 2

        #region Singularity Test 3

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty3()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T3);
            Assert.IsFalse(dsl.IsSingular());
        }

        #endregion Singularity Test 3

        #region Singularity Test 4

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty4()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T4);
            Assert.IsFalse(dsl.IsSingular());
        }

        #endregion Singularity Test 4

        #region Singularity Test 5

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty5()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Assert.IsFalse(dsl.IsSingular());
        }

        #endregion Singularity Test 5

        #region Singularity Test 10

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty10()
        {
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            Assert.IsFalse(dsl.IsSingular());
        }

        #endregion Singularity Test 10

        #region Singularity Test - Negative Case

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty()
        {
            Vector T = VectorBuilder.CreateVector(10);
            for (int i = 1; i < 10; i++)
            {
                T[i] = (double)(i + 1);
            }
            T[0] = -2.0;

            SymmetricLevinson dsl = new SymmetricLevinson(T);
            Assert.IsTrue(dsl.IsSingular());
        }

        #endregion Singularity Test - Negative Case

        #region Determinant Method Test 1

        // Test the Determinant
        [Test]
        public void DeterminantMethod1()
        {
            // calculate determinant from diagonal
            SymmetricLevinson dsl = new SymmetricLevinson(T1);

            // check results match
            Double e = System.Math.Abs((dsl.Determinant() - Det1) / Det1);
            Assert.IsTrue(e < Tolerance1);
        }

        #endregion Determinant Method Test 1

        #region Determinant Method Test 2

        // Test the Determinant
        [Test]
        public void DeterminantMethod2()
        {
            // calculate determinant from diagonal
            SymmetricLevinson dsl = new SymmetricLevinson(T2);

            // check results match
            Double e = System.Math.Abs((dsl.Determinant() - Det2) / Det2);
            Assert.IsTrue(e < Tolerance2);
        }

        #endregion Determinant Method Test 2

        #region Determinant Method Test 3

        // Test the Determinant
        [Test]
        public void DeterminantMethod3()
        {
            // calculate determinant from diagonal
            SymmetricLevinson dsl = new SymmetricLevinson(T3);

            // check results match
            Double e = System.Math.Abs((dsl.Determinant() - Det3) / Det3);
            Assert.IsTrue(e < Tolerance3);
        }

        #endregion Determinant Method Test 3

        #region Determinant Method Test 4

        // Test the Determinant
        [Test]
        public void DeterminantMethod4()
        {
            // calculate determinant from diagonal
            SymmetricLevinson dsl = new SymmetricLevinson(T4);

            // check results match
            Double e = System.Math.Abs((dsl.Determinant() - Det4) / Det4);
            Assert.IsTrue(e < Tolerance4);
        }

        #endregion Determinant Method Test 4

        #region Determinant Method Test 5

        // Test the Determinant
        [Test]
        public void DeterminantMethod5()
        {
            // calculate determinant from diagonal
            SymmetricLevinson dsl = new SymmetricLevinson(T5);

            // check results match
            Double e = System.Math.Abs((dsl.Determinant() - Det5) / Det5);
            Assert.IsTrue(e < Tolerance5);
        }

        #endregion Determinant Method Test 5

        #region Determinant Method Test 10

        // Test the Determinant
        [Test]
        public void DeterminantMethod10()
        {
            // calculate determinant from diagonal
            SymmetricLevinson dsl = new SymmetricLevinson(T10);

            // check results match
            Double e = System.Math.Abs((dsl.Determinant() - Det10) / Det10);
            Assert.IsTrue(e < Tolerance10);
        }

        #endregion Determinant Method Test 10

        #region Get Inverse 1

        // calculate inverse using Inverse member
        [Test]
        public void Inverse1()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T1);

            // check inverse
            Matrix I = dsl.Inverse();
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I1[i, j] - I[i, j]) / I1[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse 1

        #region Get Inverse 2

        // calculate inverse using Inverse member
        [Test]
        public void Inverse2()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T2);

            // check inverse
            Matrix I = dsl.Inverse();
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I2[i, j] - I[i, j]) / I2[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse 2

        #region Get Inverse 3

        // calculate inverse using Inverse member
        [Test]
        public void Inverse3()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T3);

            // check inverse
            Matrix I = dsl.Inverse();
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I3[i, j] - I[i, j]) / I3[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse 3

        #region Get Inverse 4

        // calculate inverse using Inverse member
        [Test]
        public void Inverse4()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T4);

            // check inverse
            Matrix I = dsl.Inverse();
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I4[i, j] - I[i, j]) / I4[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse 4

        #region Get Inverse 5

        // calculate inverse using Inverse member
        [Test]
        public void Inverse5()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T5);

            // check inverse
            Matrix I = dsl.Inverse();
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I5[i, j] - I[i, j]) / I5[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse 5

        #region Get Inverse 10

        // calculate inverse using Inverse member
        [Test]
        public void Inverse10()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T10);

            // check inverse
            Matrix I = dsl.Inverse();
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I10[i, j] - I[i, j]) / I10[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse 10

        #region Get Inverse Result 1

        // calculate inverse using Inverse member
        [Test]
        public void InverseResult1()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dl = new SymmetricLevinson(T1);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Inverse(I);
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I1[i, j] - I[i, j]) / I1[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse Result 1

        #region Get Inverse Result 2

        // calculate inverse using Inverse member
        [Test]
        public void InverseResult2()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dl = new SymmetricLevinson(T2);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Inverse(I);

            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I2[i, j] - I[i, j]) / I2[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }

            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse Result 2

        #region Get Inverse Result 3

        // calculate inverse using Inverse member
        [Test]
        public void InverseResult3()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dl = new SymmetricLevinson(T3);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Inverse(I);

            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I3[i, j] - I[i, j]) / I3[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse Result 3

        #region Get Inverse Result 4

        // calculate inverse using Inverse member
        [Test]
        public void InverseResult4()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dl = new SymmetricLevinson(T4);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Inverse(I);

            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I4[i, j] - I[i, j]) / I4[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse Result 4

        #region Get Inverse Result 5

        // calculate inverse using Inverse member
        [Test]
        public void InverseResult5()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dl = new SymmetricLevinson(T5);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Inverse(I);

            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I5[i, j] - I[i, j]) / I5[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse Result 5

        #region Get Inverse Result 10

        // calculate inverse using Inverse member
        [Test]
        public void InverseResult10()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dl = new SymmetricLevinson(T10);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Inverse(I);

            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I10[i, j] - I[i, j]) / I10[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse Result 10

        #region SolveVector 1

        // Test solving a linear system
        [Test]
        public void SolveVector1()
        {
            int i;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T1);
            Vector X = dsl.Solve(Y1);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((X1[i] - X[i]) / X1[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector 1

        #region SolveVector 2

        // Test solving a linear system
        [Test]
        public void SolveVector2()
        {
            int i;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T2);
            Vector X = dsl.Solve(Y2);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((X2[i] - X[i]) / X2[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector 2

        #region SolveVector 3

        // Test solving a linear system
        [Test]
        public void SolveVector3()
        {
            int i;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T3);
            Vector X = dsl.Solve(Y3);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((X3[i] - X[i]) / X3[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector 3

        #region SolveVector 4

        // Test solving a linear system
        [Test]
        public void SolveVector4()
        {
            int i;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T4);
            Vector X = dsl.Solve(Y4);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((X4[i] - X[i]) / X4[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector 4

        #region SolveVector 5

        // Test solving a linear system
        [Test]
        public void SolveVector5()
        {
            int i;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Vector X = dsl.Solve(Y5);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((X5[i] - X[i]) / X5[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector 5

        #region SolveVector 10

        // Test solving a linear system
        [Test]
        public void SolveVector10()
        {
            int i;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            Vector X = dsl.Solve(Y10);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((X10[i] - X[i]) / X10[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector 10

        #region Solve Exception Tests
        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNull()
        {
            Vector vector = null;
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            Vector X = dsl.Solve(vector);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNullResult1()
        {
            Vector vector = null;
            Vector result = VectorBuilder.CreateVector(T10.Count);
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            dsl.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNullResult2()
        {
            Vector vector = VectorBuilder.CreateVector(T10.Count);
            Vector result = null;
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            dsl.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformablel()
        {
            Vector vector = VectorBuilder.CreateVector(T10.Count + 1);
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            Vector X = dsl.Solve(vector);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformableResultNull1()
        {
            Vector vector = VectorBuilder.CreateVector(T10.Count + 1);
            Vector result = VectorBuilder.CreateVector(T10.Count);
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            dsl.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformableResultNull2()
        {
            Vector vector = VectorBuilder.CreateVector(T10.Count);
            Vector result = VectorBuilder.CreateVector(T10.Count + 1);
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            dsl.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void SolveVectorSingular1()
        {
            Vector T = VectorBuilder.CreateVector(10);
            for (int i = 1; i < 10; i++)
            {
                T[i] = (double)(i + 1);
            }
            T[0] = -2.0;

            SymmetricLevinson dsl = new SymmetricLevinson(T);

            Vector vector = VectorBuilder.CreateVector(T10.Count);
            Vector X = dsl.Solve(vector);
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void SolveVectorSingular2()
        {
            Vector vector = VectorBuilder.CreateVector(T10.Count);
            Vector result = VectorBuilder.CreateVector(T10.Count);
            Vector T = VectorBuilder.CreateVector(10);
            for (int i = 1; i < 10; i++)
            {
                T[i] = (double)(i + 1);
            }
            T[0] = -2.0;

            SymmetricLevinson dsl = new SymmetricLevinson(T);

            dsl.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNull()
        {
            Matrix matrix = null;
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            Matrix X = dsl.Solve(matrix);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNullResult1()
        {
            Matrix matrix = null;
            Matrix result = MatrixBuilder.CreateMatrix(T10.Count);
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            dsl.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNullResult2()
        {
            Matrix matrix = MatrixBuilder.CreateMatrix(T10.Count);
            Matrix result = null;
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            dsl.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformablel()
        {
            Matrix matrix = MatrixBuilder.CreateMatrix(T10.Count + 1);
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            Matrix X = dsl.Solve(matrix);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformableResultNull1()
        {
            Matrix matrix = MatrixBuilder.CreateMatrix(T10.Count + 1);
            Matrix result = MatrixBuilder.CreateMatrix(T10.Count);
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            dsl.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformableResultNull2()
        {
            Matrix matrix = MatrixBuilder.CreateMatrix(T10.Count);
            Matrix result = MatrixBuilder.CreateMatrix(T10.Count + 1);
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            dsl.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformableResultNull3()
        {
            Matrix matrix = MatrixBuilder.CreateMatrix(T10.Count, T10.Count);
            Matrix result = MatrixBuilder.CreateMatrix(T10.Count, T10.Count + 1);
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            dsl.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void SolveMatrixSingular1()
        {
            Vector T = VectorBuilder.CreateVector(10);
            for (int i = 1; i < 10; i++)
            {
                T[i] = (double)(i + 1);
            }
            T[0] = -2.0;

            SymmetricLevinson dsl = new SymmetricLevinson(T);

            Matrix matrix = MatrixBuilder.CreateMatrix(T10.Count);
            Matrix X = dsl.Solve(matrix);
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void SolveMatrixSingular2()
        {
            Matrix matrix = MatrixBuilder.CreateMatrix(T10.Count);
            Matrix result = MatrixBuilder.CreateMatrix(T10.Count);
            Vector T = VectorBuilder.CreateVector(10);
            for (int i = 1; i < 10; i++)
            {
                T[i] = (double)(i + 1);
            }
            T[0] = -2.0;

            SymmetricLevinson dsl = new SymmetricLevinson(T);

            dsl.Solve(matrix, result);
        }

        #endregion

        #region SolveVector Result 1

        // Test solving a linear system
        [Test]
        public void SolveVectorResult1()
        {
            int i;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T1);
            Vector X = VectorBuilder.CreateVector(T1.Count);
            dsl.Solve(Y1, X);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((X1[i] - X[i]) / X1[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector Result 1

        #region SolveVector Result 2

        // Test solving a linear system
        [Test]
        public void SolveVectorResult2()
        {
            int i;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T2);
            Vector X = VectorBuilder.CreateVector(T2.Count);
            dsl.Solve(Y2, X);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((X2[i] - X[i]) / X2[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector Result 2

        #region SolveVector Result 3

        // Test solving a linear system
        [Test]
        public void SolveVectorResult3()
        {
            int i;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T3);
            Vector X = VectorBuilder.CreateVector(T3.Count);
            dsl.Solve(Y3, X);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((X3[i] - X[i]) / X3[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector Result 3

        #region SolveVector Result 4

        // Test solving a linear system
        [Test]
        public void SolveVectorResult4()
        {
            int i;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T4);
            Vector X = VectorBuilder.CreateVector(T4.Count);
            dsl.Solve(Y4, X);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((X4[i] - X[i]) / X4[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector Result 4

        #region SolveVector Result 5

        // Test solving a linear system
        [Test]
        public void SolveVectorResult5()
        {
            int i;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T5);
            Vector X = VectorBuilder.CreateVector(T5.Count);
            dsl.Solve(Y5, X);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((X5[i] - X[i]) / X5[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector Result 5

        #region SolveVector Result 10

        // Test solving a linear system
        [Test]
        public void SolveVectorResult10()
        {
            int i;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T10);
            Vector X = VectorBuilder.CreateVector(T10.Count);
            dsl.Solve(Y10, X);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                e = System.Math.Abs((X10[i] - X[i]) / X10[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector Result 10

        #region Solve Matrix 1

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrix1()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T1);

            // check inverse
            Matrix I = dsl.Solve(MatrixBuilder.CreateIdentityMatrix(1));
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I1[i, j] - I[i, j]) / I1[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix 1

        #region Solve Matrix 2

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrix2()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T2);

            // check inverse
            Matrix I = dsl.Solve(MatrixBuilder.CreateIdentityMatrix(2));
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I2[i, j] - I[i, j]) / I2[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix 2

        #region Solve Matrix 3

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrix3()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T3);

            // check inverse
            Matrix I = dsl.Solve(MatrixBuilder.CreateIdentityMatrix(3));
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I3[i, j] - I[i, j]) / I3[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix 3

        #region Solve Matrix 4

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrix4()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T4);

            // check inverse
            Matrix I = dsl.Solve(MatrixBuilder.CreateIdentityMatrix(4));
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I4[i, j] - I[i, j]) / I4[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix 4

        #region Solve Matrix 5

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrix5()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T5);

            // check inverse
            Matrix I = dsl.Solve(MatrixBuilder.CreateIdentityMatrix(5));
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I5[i, j] - I[i, j]) / I5[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix 5

        #region Solve Matrix 10

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrix10()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T10);

            // check inverse
            Matrix I = dsl.Solve(MatrixBuilder.CreateIdentityMatrix(10));
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I10[i, j] - I[i, j]) / I10[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix 10

        #region Solve Matrix Result 1

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrixResult1()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T1);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(T1.Count);
            dsl.Solve(MatrixBuilder.CreateIdentityMatrix(1), I);
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I1[i, j] - I[i, j]) / I1[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix Result 1

        #region Solve Matrix Result 2

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrixResult2()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T2);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(T2.Count);
            dsl.Solve(MatrixBuilder.CreateIdentityMatrix(2), I);
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I2[i, j] - I[i, j]) / I2[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix Result 2

        #region Solve Matrix Result 3

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrixResult3()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T3);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(T3.Count);
            dsl.Solve(MatrixBuilder.CreateIdentityMatrix(3), I);
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I3[i, j] - I[i, j]) / I3[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix Result 3

        #region Solve Matrix Result 4

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrixResult4()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T4);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(T4.Count);
            dsl.Solve(MatrixBuilder.CreateIdentityMatrix(4), I);
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I4[i, j] - I[i, j]) / I4[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix Result 4

        #region Solve Matrix Result 5

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrixResult5()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T5);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(T5.Count);
            dsl.Solve(MatrixBuilder.CreateIdentityMatrix(5), I);
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I5[i, j] - I[i, j]) / I5[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix Result 5

        #region Solve Matrix Result 10

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrixResult10()
        {
            int i, j;
            double e, me;
            SymmetricLevinson dsl = new SymmetricLevinson(T10);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(T10.Count);
            dsl.Solve(MatrixBuilder.CreateIdentityMatrix(10), I);
            me = 0.0;
            for (i = 0; i < dsl.Order; i++)
            {
                for (j = 0; j < dsl.Order; j++)
                {
                    e = System.Math.Abs((I10[i, j] - I[i, j]) / I10[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix Result 10
    }
}
